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Abstract. In this paper we present a novel fast method to solve Poisson equation in an arbitrary 
two dimensional region with Neumann boundary condition. The basic idea is to solve the original 
Poisson problem by a two-step procedure: the first one finds the electric displacement field D and 
the second one involves the solution of potential <j>. The first step exploits loop-tree decomposition 
technique that has been applied widely in integral equations within the computational electromag- 
netics (CEM) community. We expand the electric displacement field in terms of a tree basis. Then, 
coefficients of the tree basis can be found by the fast tree solver in O(N) operations. Such ob- 
tained solution, however, fails to expand the exact field because the tree basis is not completely curl 
free. Despite of this, the accurate field could be retrieved by carrying out a procedure of divergence 
free field removal. Subsequently, potential distribution <j> can be found rapidly at the second stage 
with another fast approach of O(N) complexity. As a result, the method dramatically reduces so- 
lution time comparing with traditional FEM with iterative method. Numerical examples including 
electrostatic simulations are presented to demonstrate the efficiency of the proposed method. 
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1. Introduction. There are a variety of physical situations and engineering 
problems described by elliptic partial differential equations (PDEs) such as the Pois- 
son equation: 

V 2 u(r) = /(r). 

Examples of this equation are encountered in low-frequency dielectric or conductivity 
problems PQPJ. This equation is often solved in micro and nanoelectronic device 
physics: it is also found in electronic transport and electrochemistry in terms of the 
Poisson-Boltzmann equation [3] [4] . Hence, finding a robust and efficient way to solve 
it has attracted considerable interests in various fields of research. 

Over the past few decades, several kinds of fast methods for solving Poisson 
equation have been proposed. One popular fast Poisson solver is based on Fourier 
analysis and accelerated by FFT [2]. However, this method has generally been limited 
to regular geometries, such as rectangular regions, 2D polar and spherical geometries 
[5], and spherical shells [B]. 

Multigrid methods are generally accepted as among the fastest numerical methods 
[7-9] . These methods take advantage of fine mesh and coarse mesh. Multigrid methods 
can be used in irregular domains and also be extended to other partial differential 
equations. But it is difficult to implement multigrid methods in a robust fashion 
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Fig. 2.1. Schema of the electrostatic problem 



because they demand a hierarchy of grids of different density, which are not convenient 
in many real world problems. 

Another algorithm known as the fast multipole method (FMM) [10-12] has been 
developed to solve Poisson problems with 0(N log N) complexity. This method is 
applied to integral equation derived via the Green's function rather than differential 
methods where Poisson equation is discritized directly. FMM solvers are particularly 
well suited for solving irregular shape problems. 

In this paper, we will solve Poisson equation with Neumann boundary condition, 
which is often encountered in electrostatic problems, through a newly proposed fast 
method. Instead of discretizing Poisson equation directly, we solve it in two sequential 
steps: The first step aims to find the displacement electric field D from V • D = p, 
where p is the charge distribution; then, the second one obtains potential <j) from 
V0 = — D/e, in which e stands for the permittivity of the medium. 

The remainder of the paper is organized as follows. In Section [2] we outline the 
problem arising from electrostatics that we seek to solve. In Section [3j details of the 
proposed method will be described, in which two sequential steps of this method will 
be presented by two subsections. Next, in Section[4j we demonstrate and validate the 
efficiency by applying several numerical examples. Finally, the conclusions are drawn 
in Section [5] 

2. Problem Description. The general Poisson equation for heterogeneous me- 
dia is expressed as: 



where </>(r) is the potential and p(r) describes the charge density, e r (r) and eo are 
relative permittivity and free space permittivity, respectively. 

In this paper, we focus on this kind of Poisson problems, as illustrated in Fig. l2.I[ 
where the charge density p{x, y) distributes within a two dimensional domain SI and 
the relative permittivity e r is a constant. Furthermore, the Neumann boundary con- 
dition is imposed on T. 

Thus the equation to be solved is 



(2.1) 



V • e r (r)V0(r) 



P(r) 



for (a:, y)e(! 



e r eo 



(2.2) 



dcf)(x,y) 



for (x, y) 6 r. 
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3. Solution method. The proposed new method solves (|2.2j) by two steps. 
More specifically, we first solve this equation 

(3.1) W-T> = p(x,y) 

where the electric displacement field D is the product of permittivity e and electric 
field E, namely 

(3.2) D = eE 
and 

(3.3) e = e r e . 

Once we have the electric displacement field, the potential distribution <j) can be found 
by solving 

(3.4) V4>(x,y) = --, 

e 

which corresponds to the second step. 

Although two steps are required for this algorithm, fast methods can be applied 
to both steps by virtue of the tree basis. 

3.1. Solution of the Electric Field. In order to find the electric field distri- 
bution, the problem at this stage is represented by 

V • D = p(x, y) for (x, y) G CI 

(3.5) n.D = ^^ for (x,y) eT 

e r eo 

and 

(3.6) VxE = 0. 

For the homogeneous medium, we can further infer that the electric displacement 
field is curl free as well (V x D = 0), since the relative permittivity e r is constant. 

3.1.1. Construct Matrix System by the Tree Basis. In computational elec- 
tromagnetics community, the low frequency problem has attracted intensive research 
for the last decade. When the frequency is low, the current J naturally decomposes 
into a solenoidal (divergence free) part and an irrotational (curl free) part which 
is known as the Hclmholtz decomposition [13]. Moreover, these two parts are not 
balanced when the frequency becomes low. Hence, there is a severe numerical prob- 
lem when solving integral equations in which RWG bases are normally used. One 
remarkable remedy to this low frequency breakdown is the well known loop-tree de- 
composition [14-18]. The RWG basis set is decomposed into the loop basis, which 
has zero divergence, and the tree basis, which has non-zero divergence. This is a 
quasi-Helmholtz decomposition because the tree basis is not curl free. 

Borrowing the idea of loop-tree decomposition, we can use a tree basis to expand 
the irrotational electric displacement field as 

Nt 

(3.7) T> tree = J2 t i T i( x ,y) 

i=l 
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where Tj is the i-th tree basis whose coefficient is tj, and N t is the total number of 
tree basis functions. It is equal to the number of patches minus one for a triangular 
mesh. We then use the pulse basis to expand the charge density p, namely 

(3.8) p(x, y) = ^CjPj (x,y) 

i=l 

where N p stands for the patch number. The pulse function is defined as 

p , . I 1 if (x, y) G i-th patch 
%\ x i V) - | q otherwise 

Substituting (|3.7p and (|3.8[) into p. 51) and testing it with the same set of pulse 
functions as in the Galerkin's method, we have a matrix equation 

(3.9) K I t = V p 
where 

(3.10) [K\..= J J P i (x,y)V-T j (x,y)dxdy 



(3.H) ^ Vp ^' = y / Pi ( x > y } p ( x > y "> dxdy = Ci 

3.1.2. Fast Tree Solver. As described in [13] [16], (O can be solved in 0{N) 
operations. We outline the basic idea as follows. In the first place, a tree is found to 
span the surface where Poisson equation is solved. This tree corresponds to a set of 
RWG basis functions complementary to the loop basis. Fig. 13.11 shows a typical tree 
structure. Each node of the tree represents the center of a patch, which is related to 
the triangle patch of the RWG basis. Each line between two nodes can be associated 
with the current that flows between two patches. Furthermore, the patch unknown 
associated with a tip of the dendritic branch is only related to one current unknown. 
Therefore, starting from the branch tips, the unknowns can be solved for recursively 
until a junction is reached. Noticing the fact that a junction node cannot connect 
more than three neighboring nodes for the case of RWG functions, we can solve for 
the unknowns of junctions when the unknowns on two associated open branches have 
been solved. Hence, all unknowns can be solved for in O(N) operations. 

3.1.3. Divergence Free Field Removal. It is well known that the tree ba- 
sis is not curl free. This property makes the D tree from (|3.9[) inaccurate because 
some solenoidal components contaminate it. Moreover, the tree basis is non-unique. 
Theoretically, the ~D tr ee is a vector living in the space spanned by the RWG basis, 
which satisfies Eq. (|3.5j) . From this point of view, the desired D can be obtained by 
projecting Dtree onto the curl free space. However, it is formidable to find a basis set 
which is completely curl free in the RWG space. In [19] , SVD method was used for 
finding the solenoidal basis set. But SVD is expensive in terms of both memory and 
computing time requirements. 

Another way to circumvent this obstacle is to remove the complementary solenoidal 
(divergence free) part that is orthogonal to the curl free space. This is practical be- 
cause the divergence free part, unlike its curl free counterpart, could be expanded by 
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Fig. 3.1. Schematic of a general tree 



the loop basis that is living on the space spanned by RWG functions. It is noted that 
the loop basis has divergence free property; hence, it can expand the divergence free 
space completely. Therefore, we can first use a loop basis to project D tree onto the 
divergence free space. Once the projection has been done, the pure curl free part of 
Dtree can be obtained by subtracting the divergence free components. 
Mathematically, Dtree can be expressed by 



N 



Ni 



(3.12) 



Dtree = ^V^^T,?/) + ^ Ijlij (x, y) 



where the first term on right side corresponds to pure curl free part of D tree , while 
the second one associates the divergence free part. Moreover, refers to RWG basis 
functions whose total number is TV, and Lj denotes the loop basis functions with total 
number of JV/. 

By using the loop basis to test the above equation, Eq. (|3.12p leads to a matrix 
system 



(3.13) 
where 



Gfli= V 



(3.14) t G 4=/ Ju(x,y)-L J (x,y)dxdy 
is the Gram matrix of loop basis, and 

(3.15) [Vali = J J Li(x,y) ■D tree (x,y)dxdy 

is the projection value of D tree onto the loop space. In this equation, we use the 
relation 



(3.16) 



n 



Li(x, y) ■ ^2 a &i x ' v) & x &y = 



i=l 



() 
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because the fact that loop basis is orthogonal to bases of curl free space. 

The Gram matrix of loop bases, G/, is highly sparse. Moreover, it is a symmetric, 
positive definite and diagonal dominant matrix. Therefore, common iterators, such 
as CG, Bi-CGSTAB [3D] and GMRES HQ [22], work efficiently for solving Eq. (|3TT3"]) . 

Consequently, the electric displacement field desired can be retrieved by 

(3.17) D = D tree -^Z i L i 

i=i 

3.2. Solution of Potential. In the second stage, we need to solve the following 
equation to find potential 

(3.18) V<j>(x,y) = --. 

e 

This can be achieved by using the same method of the above section because the del 
operator (V) is the transpose of the divergence operator (V-). Hence, we expand the 
potential <fi in terms of pulse functions, that is, 

(3.19) <f>{x,y)=Y^v i P i (x,y) 

i=l 

With the similar method described in Section [3. 11 testing this equation using a set of 
tree basis, we obtain 

(3.20) K T I = V 
where 

(3.21) [V^ = - / [Tfay). dxdy 

and = Vi which is defined in fl3.19p . The element K is originally defined as 



— T 

K 



J J Ti(x,y)-VPj(x,y)dxdy, 



which can be deduced using integration by parts as 



(3.22) [K j = J J P,(x,y)\7 ■ T t (x, y) dxdy = [K] .. 

Consequently, the resulting matrix is just a transpose of the matrix [K] from 
Section [3~T1 Hence, (I3.20P can be solved with fast tree solver in O(N) operations. 

4. Numerical Examples. In this section, two numerical examples are shown 
to validate the efficiency of the proposed method. All examples have been calculated 
on a standard computer with 2.66 GHz CPU, 4 GB memory and Windows operating 
system. 
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(a) Simulation potential distribution 



(b) Exact potential distribution 



Fig. 4.1. Potential distribution\(a)\ Simulation result\(bj\ Exact potential distribution 



4.1. Simple Neumann Problem. In order to validate the correctness of this 
algorithm, we solve a simple 2-D Poisson equation as the first example. The Poisson 
equation in this case is 



(4.1) 



— -z + tt^ — — 7rcos(7ra;) - 7rcos(7rj/J 
ox z oy z 



where (x,y) 6 f! = [0, 1] X [0, 1], with the Neumann boundary condition 



(4.2) 



dx 



dy 







[x,y)eT. 



It is well known that Neumann problems is uniquely solvable but up to a con- 
stant; hence, we have to impose a reference potential so that the uniqueness can be 
guaranteed. By setting the potential as 2/tt at the point (0,0), the aforementioned 
problem has an analytical solution 



(4.3) 



4>(x, y) = [cos(7rx) + cos(7ry)] /it. 



By utilizing the proposed method, potential distribution can be found as shown 



in Fig. 14.11 where Fig. 4.1(a) shows the potential distribution result while Fig. 4.1(b) 
gives the exact solution as a reference. Furthermore, in Fig. 14.21 we compared the 
value from the proposed approach with the analytic results. It is obvious that the 
potential distribution from this new method agrees with analytical one well. 

4.2. Line Source Problem. As the second example, a Poisson problem involv- 
ing two line charge sources was simulated within a polygonal region. One line source 
was at point (—0.475, —0.329) with a positive unit charge while the other one was at 
point (0.494, —0.459) with a negative unit charge. The resultant potential distribution 
is shown in Fig. 14.31 In this case, the potential value of the up-left corner (—2.0, 1.0) 
was set to be zero as a reference potential. The result has been validated by FEM 
method. 

4.3. Computational Complexity Analysis. To analyze the computational 
complexity, we have computed the problem of Section ^. ll with different mesh densities. 
The solution time of this algorithm consists of three parts: the first fast tree solution 
time for (|3.9|) . the divergence free field removal time and the second fast tree solution 
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Fig. 4.3. Potential distribution calculated by the proposed method for line sources 



time for (|3.20p . As for the fast tree solution time, Zhao and Chew have proved that 
it is of O(N) complexity [TB]. Hence, the bottleneck comes from the procedure of 
divergence free field removal, which amounts to the iterative solution time of (|3.13p . 
The CPU time of this part therefore depends on iterative solver type and the required 
accuracy level. 

In divergence free field removal part, the Gram matrix of loop basis has the form 
(4.4) [G,] tj = J J Li{x,y) ■ L 3 (x,y) dxdy. 

As similar to the stiffness matrix of FEM, this matrix is a symmetric, positive definite 
system. For this kind of matrix systems, the number of iterations is proportional to 
the square-root of the condition number. 

Fig. 14.41 shows the computational complexity comparison between this fast Pois- 
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Fig. 4.4. Complexity comparison with same stopping criterion 



son solver and FEM. Both conjugated gradient (CG) and GMRES are adopted with 
same stopping criterion (0.01). For GMRES, the restart parameter is 60. As can be 
seen from this plot, the total solution time of our proposed new method is evidently 
less than the one of traditional FEM, whichever CG or GMRES is adopted. In our nu- 
merical experiment, the complexity of this new method approaches linear complexity. 
Moreover, GMRES methods work better than CG method in our simulations. 

5. Conclusions. We have proposed and developed a new two dimensional fast 
Poisson solver for electrostatic problems with the Neumann boundary condition based 
on the loop-tree decomposition technique. This method could solve Poisson equation 
rapidly with the fast tree solver. Its computational consumption is evidently less than 
the one of traditional FEM. This method promises to be a novel fast Poisson solver 
that can solve general Poisson equation for electrostatics as well as other related fields. 

Acknowledgments. The authors would like to thank Dr. Shcng Sun, Dr. Min 
Tang and Jun Huang for their helpful discussions and feedback from Leslie Greengard. 



REFERENCES 



[1] J. D. Jackson, Classical Electrodynamics Third Edition, Third ed., Wiley, 1998. 

[2] S. N. BARKAS, An introduction to fast poisson solvers, 

http:/ /people. freebsd.org/~snb/school/fastpoisson.pdf 2005. 
[3] F. FOGOLARI, A. Brigo, AND H. Molinari, The Poisson-Boltzmann equation for biomolecular 

electrostatics: a tool for structural biology, Journal of Molecular Recognition, 15 (2002), 

pp. 377-392. 

[4] A. Adelmann, P. Arbenz, and Y. Ineichen, A fast parallel poisson solver on irregular do- 
mains applied to beam dynamics simulations, Journal of Computational Physics, 229(2010), 
pp. 4554 - 4566. 

[5] M. Lai and W. Wang, Fast direct solvers for poisson equation on 2D polar and spherical 
geometries, Numerical Methods for Partial Differential Equations, 18(2002), pp. 56—68. 

[6] Y.-L. Huang, J.-G. Liu, and W.-C. Wang, An FFT based fast poisson solver on spherical 
shells, Communications in Computational Physics, 9(2011), pp. 649—667. 

[7] U. Trottenberg, C. W. Oosterlee and A. Schuller, Multigrid, Academic Press, 2001. 

[8] S. R. Fulton, P. E. Ciesielski, and W. H. Schubert, Multigrid methods for elliptic problems: 
A review, Monthly Weather Review, 14(1986), pp. 943-959. 



10 A NOVEL FAST SOLVER FOR POISSON EQUATION 

[9] L. Briggs, V. E. Henson, and S. F. McCormick, A Multigrid Tutorial, Philadelphia: SIAM, 
2000. 

[10] A. McKenney, L. Greengard, and A. MAYO, A fast poisson solver for complex geometries, 
Journal of Computational Physics, 118(1995), pp. 348 - 355. 

[11] F. Ethridge and L. Greengard, A new fast-multipole accelerated poisson solver in two di- 
mensions, SIAM J. on Scientific Computing, 23(2001), pp. 741-760. 

[12] M. H. Langston, L. Greengard, and D. Zorin, A free-space adaptive fmm-based pde solver 
in three dimensions, Communications in Applied Mathematics and Computational Science, 
6(2011), pp. 79-122. 

[13] W. C. Chew, M. S. Tong, and B. Hu, Integral Equations Methods for Electromagnetic and 
Elastic Waves, Morgan & Claypool, 2008. 

[14] D. R. Wilton and A. W. Glisson, On improving the electric field integral equation at low 
frequencies, 1981 Spring URSI Radio Science Meeting Digest, Los Angeles, CA, 1981, p. 24. 

[15] J. Mautz and R. Harrington, An E-field solution for a conducting surface small or compa- 
rable to the wavelength, Antennas and Propagation, IEEE Transactions on, 32(1984), pp. 
330 - 339. 

[16] J.-S. Zhao and W. C. Chew, Integral equation solution of maxwell's equations from zero 

frequency to microwave frequencies, Antennas and Propagation, IEEE Transactions on, 

48(2000), pp. 1635 -1645. 
[17] W. Wu, A. W. Glisson, and D. Kajfez, A comparison of two low-frequency formulations 

for the electric field integral equation, Tenth Ann. Rev. Prog. Appl. Comput. Elcctromag., 

2(1994), pp. 484-491. 

[18] M. Burton and S. Kashyap, A study of a recent, moment-method algorithm that is accurate 

to very low frequencies, Appl. Comput. Elcctromagn. Soc. J., 10(1995), pp. 58—68. 
[19] F. Vipiana, P. Pirinoli, and G. Vecchi, A multiresolution method of moments for triangular 

meshes, Antennas and Propagation, IEEE Transactions on, 53(2005), pp. 2247 - 2258. 
[20] H. A. VAN der VORST, Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the 

solution of nonsymmetric linear systems, SIAM J. on Scientific Computing, 13(1992), pp. 

631-644. 

[21] Y. Saad, Iterative Methods for Sparse Linear Systems, Second ed., Society for Industrial and 

Applied Mathematics, 2003. 
[22] Y. Saad and M. Schultz, GMRES: A generalized minimal residue algorithm for solving 

nonsymmetric linear systems, SIAM J. Sci. Stat. Comput., 7(1986), pp. 856-869. 



